#suppressWarnings(library(bitops))
#suppressWarnings(library(splines))
#suppressWarnings(library(gtools))
#suppressWarnings(library(gdata))
#suppressWarnings(library(caTools))
#suppressWarnings(library(grid))
#suppressWarnings(library(survival))
#suppressWarnings(library(gplots))
#suppressWarnings(library(mvtnorm))
#suppressWarnings(library(multcomp))
#suppressWarnings(library(car))
suppressWarnings(library(R2HTML))
#suppressWarnings(library(lattice))
suppressWarnings(library(nlme))

# retrieve args
Args <- commandArgs(TRUE)

#Read in data
statdata <- read.csv(Args[3], header=TRUE, sep=",")

#Copy Args


responseTransform <- Args[5]
covariateTransform <- Args[6]
treatments <- Args[7]
otherFactors <- Args[8]
covariate <- Args[9]
sigLevel <- Args[10]
randomEffect1 <- Args[11]
randomEffect2 <- Args[12]
randomEffect3 <- Args[13]
randomEffect4 <- Args[14]
repList1 <- Args[15]
repList2 <- Args[16]
repList3 <- Args[17]
repList4 <- Args[18]


#STB Edits so that covariate can be non-null but treatments are null -  17Oct2010
if (treatments=="NULL" && covariate == "NULL")
	{
	nonetreat<-eval(parse(text = paste(Args[4], "1")))
	print(nonetreat)
	model <- as.formula(nonetreat)
	model2 <- nonetreat
	} else
	{
	model <- as.formula(Args[4])
	model2 <- Args[4]
	} 




#if (treatments!="NULL")
#	{
#	model <- as.formula(Args[4])
#	model2 <- Args[4]
#	} else 
#	{
#	nonetreat<-eval(parse(text = paste(Args[4], "1")))
#	print(nonetreat)
#	model <- as.formula(nonetreat)
#	model2 <- nonetreat
#	}
#

#model
#responseTransform
#covariateTransform
#treatments
#otherFactors
#randomEffect1
#randomEffect2
#randomEffect3
#randomEffect4
#repList1
#repList2
#repList3
#repList4

# Making sure random factors are factor

if (randomEffect1 != "NULL")
{
	statdata$rEffect1 <-as.factor(eval(parse(text = paste("statdata$", randomEffect1))))
}
if (randomEffect2 != "NULL")
{
	statdata$rEffect2  <-as.factor(eval(parse(text = paste("statdata$", randomEffect2))))
}
if (randomEffect3 != "NULL")
{
	statdata$rEffect3  <-as.factor(eval(parse(text = paste("statdata$", randomEffect3))))
}
if (randomEffect4 != "NULL")
{
	statdata$rEffect4  <-as.factor(eval(parse(text = paste("statdata$", randomEffect4))))
}


#parameters sorted in the R code

#covariate<-"covariate1"

#statdata$alltreat<-statdata$Treatment
siglevel<-as.numeric(sigLevel)

statdata$interactionEffect<-NULL


#Calculating the highest order interaction

if (treatments!="NULL")
	{
	treatments1 <- unlist(strsplit(treatments,","))
	lentr<-length(treatments1)
	firsttreat<-treatments1[1]
	statdata$alltreat <-eval(parse(text = paste("statdata$", firsttreat)))

	if (lentr >= 2)
		{
			for (i in 2:lentr)
			{
			secondtreat<-treatments1[i]
			statdata$alltreat2 <-eval(parse(text = paste("statdata$", secondtreat)))
			alltreat<-paste(statdata$alltreat, " ",statdata$alltreat2)
			statdata$alltreat<-NULL
			alltreat<-as.factor(alltreat)
			statdata<-cbind(statdata,alltreat)
			}
		}
	} else {
		statdata$alltreat<- 1
		}



#Define random model
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL")
{
	randmodel<-as.formula(~1|rEffect1)
}
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL")
{
	randmodel<-as.formula(~1|rEffect1/rEffect2)
}
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL")
{
	randmodel<-as.formula(~1|rEffect1/rEffect2/rEffect3)
}
#randmodel<-as.formula(~1|Animal/Slide)
#randmodel<-as.formula(~1|Animal/Slide/ROI)
#randmodel<-as.formula(~1|Random1/Random2/Random3/Random4)


#Simons power analysis function
sipowerttest<-function(delta,sd,n,sig.level)
{
	NCP<-delta/(sqrt(sd*2/n))
	crit1<-qt((1-sig.level/2),df=(2*n-2))
	crit2<-qt((1-(1-sig.level/2)),df=(2*n-2))
	suppressWarnings(power<-1-pt(crit1,(2*n-2),NCP)+pt(crit2,(2*n-2),NCP))
}



#Parameters set up
resp <- unlist(strsplit(Args[4],"~"))[1] #get the response variable from the main model

graphTitle0 = "Power curves for original design"
graphTitle1 = "Power curves for predicted design"


#Setup the html file and associated css file
htmlFile <- sub(".csv", ".html", Args[3]); #determine the file name of the html file
.HTML.file = htmlFile
cssFile <- "r2html.css"
cssFile <- paste("'",cssFile,"'", sep="") #need to enclose in quotes when path has spaces in it
HTMLCSS(CSSfile = cssFile)

#Output HTML header
HTML.title("<bf>SilveR Nested Design Analysis", HR=1, align="left")




# Checking the conditions

#Only run analysis with two or more random effects
if(randomEffect1=="NULL")
{
messag1<-c("Unfortunately you need two nested random effects in order to use this module.")
print(messag1)
HTML.title("</bf> ", HR=2, align="left")
HTML.title(messag1, HR=0, align="left")
quit()
} 

if(randomEffect2=="NULL")
{
messag1<-c("Unfortunately you need two nested random effects in order to use this module.")
print(messag1)
HTML.title("</bf> ", HR=2, align="left")
HTML.title(messag1, HR=0, align="left")
quit()
} 


if(repList1=="NULL" && repList2 =="NULL" && repList3 =="NULL" && repList4 =="NULL" )
{
messag1<-c("You need to define at least one list of random effect replications to proceed with this analysis.")
print(messag1)
HTML.title("</bf> ", HR=2, align="left")
HTML.title(messag1, HR=0, align="left")
quit()
} 



#Check lowest order Random effect is defines observational unit

len<-dim(statdata)[1]
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL")
{
	len2<-length(unique(eval(parse(text = paste("statdata$", randomEffect2)))))
}
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL")
{
	len2<-length(unique(eval(parse(text = paste("statdata$", randomEffect3)))))
}
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL")
{
	len2<-length(unique(eval(parse(text = paste("statdata$", randomEffect4)))))
}
if(len2<len)
{
messag2<-c("The levels of the lowest order random effect do not define the observations units (rows of the dataset). Please include an extra random effect that does.")
print(messag2)
HTML.title("</bf> ", HR=2, align="left")
HTML.title(messag2, HR=0, align="left")
quit()
}





#Warning
title<-c("Warning")
HTML.title(title, HR=2, align="left")

HTML.title("</bf> ", HR=2, align="left")
HTML.title("Warning: This module is currently under construction, care should be taken when considering the results. The results have not been verified.", HR=0, align="left")


#if (changesValue!="NULL")
#{
#message3<-c("The Expected Changes options are not currently implemented. Any user selections will be ignored.")
#	HTML.title("</bf> ", HR=2, align="left")
#	HTML.title(message3, HR=0, align="left")
#	print(message3)
#}


#Response

HTMLbr()
title<-c("Response")
if(covariate != "NULL")
{
	title<-paste(title, " and covariate", sep="")
}

HTML.title(title, HR=2, align="left")
add<-paste(c("The  "), resp, sep="")
add<-paste(add, " response is currently being analysed by the Nested Design Analysis module", sep="")

if(covariate != "NULL")
{
	add<-paste(add, c(", with  "), sep="")
	add<-paste(add, covariate, sep="")
	add<-paste(add, " fitted as a covariate.", sep="")
} else {
	add<-paste(add, ".", sep="")
}

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")

if (responseTransform != "None" | covariateTransform != "None")
{
	HTML.title("<bf> ", HR=2, align="left")
}

if (responseTransform != "None")
{
	add2<-paste(c("The response has been "), responseTransform, sep="")
	add2<-paste(add2, " transformed prior to analysis.", sep="")
	HTML.title(add2, HR=0, align="left")
}

if (covariateTransform != "None")
{
	add3<-paste(c("The covariate has been "), covariateTransform, sep="")
	add3<-paste(add3, " transformed prior to analysis.", sep="")
	HTML.title(add3, HR=0, align="left")
}




#Method

HTMLbr()
title<-c("Methodology")
HTML.title(title, HR=2, align="left")

add<-c("This module uses the estimated variance components from the original dataset to predict the hypothetical statistical power than can be achieved by varying the replication of the random effects in the experimental design, see Snedecor and Cochran (1989, p239).")
HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")



#Model fit
modelfit<-lme(model, data=statdata, random=randmodel) 
print(anova(modelfit))
#Variance components

#Defined random effects
cov<-VarCorr(modelfit) 

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL")
{
	cov11<-as.numeric(cov[1])
	cov12<-as.numeric(cov[2])
	covs<-c(cov11, cov12)
	reffects<-c(randomEffect1, randomEffect2)
}
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL")
{
	cov21<-as.numeric(cov[2])
	cov22<-as.numeric(cov[4])
	cov23<-as.numeric(cov[5])
	covs<-c(cov21, cov22, cov23)
	reffects<-c(randomEffect1, randomEffect2, randomEffect3)
}
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL")
{
	cov31<-as.numeric(cov[2])
	cov32<-as.numeric(cov[4])
	cov33<-as.numeric(cov[6])
	cov34<-as.numeric(cov[7])
	covs<-c(cov31, cov32, cov33, cov34)
	reffects<-c(randomEffect1, randomEffect2, randomEffect3, randomEffect4)
}



#Variance components table


HTMLbr()
title<-c("Table of estimated variance components")
HTML.title(title, HR=2, align="left")

table<-matrix(nrow=1, ncol=length(reffects))
for ( i in 1:length(reffects))
{
	table[1,i]=format(round(covs[i], 4), nsmall=4, scientific=FALSE)
}
vars<-data.frame(table)
colnames(vars)<-reffects
rownames(vars)<-c("Variance Components")

HTML(vars, classfirstline="second", align="left")






#Default Replication

#Random1 replication
length1<-length(unique(levels(as.factor(statdata$alltreat))))
vectorN1 <-c(1:length1)
for (i in 1:length1)
{
	sub<-subset(statdata, statdata$alltreat == unique(levels(as.factor(statdata$alltreat)))[i])
	sub1<-data.frame(sub)
	vectorN1[i]=dim(sub1)[1]
}
repz1<-as.numeric(mean(vectorN1))


#Random2 replication
if (randomEffect1 != "NULL")
{
	length2<-length(unique(levels(as.factor(statdata$rEffect1))))
	vectorN2 <-c(1:length2)
	for (i in 1:length2)
	{
		sub<-subset(statdata, statdata$rEffect1 == unique(levels(as.factor(statdata$rEffect1)))[i])
		sub2<-data.frame(sub)
		vectorN2[i]=dim(sub2)[1]
	}
	repz2<-as.numeric(mean(vectorN2))
}

#Random3 replication
if (randomEffect2 != "NULL")
{
	length3<-length(unique(levels(as.factor(statdata$rEffect2))))
	vectorN3 <-c(1:length3)
	for (i in 1:length3)
	{
		sub<-subset(statdata, statdata$rEffect2 == unique(levels(as.factor(statdata$rEffect2)))[i])
		sub3<-data.frame(sub)
		vectorN3[i]=dim(sub3)[1]
	}
	repz3<-as.numeric(mean(vectorN3))
}

#Random4 replication
if (randomEffect3 != "NULL")
{
	length4<-length(unique(levels(as.factor(statdata$rEffect3))))
	vectorN4 <-c(1:length4)
	for (i in 1:length4)
	{
		sub<-subset(statdata, statdata$rEffect3 == unique(levels(as.factor(statdata$rEffect3)))[i])
		sub4<-data.frame(sub)
		vectorN4[i]=dim(sub4)[1]
	}
	repz4<-as.numeric(mean(vectorN4))
}



#Calculating replication for 4 random effects
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL")
{
	rep1<-as.numeric(format(round(repz1/repz2, 0), nsmall=0, scientific=FALSE))
	rep2<-as.numeric(format(round(repz2/repz3, 0), nsmall=0, scientific=FALSE))
	rep3<-as.numeric(format(round(repz3/repz4, 0), nsmall=0, scientific=FALSE))
	rep0<-repz4
	reps3<-c(rep1, rep2, rep3, rep0)
	reps<-reps3
}

#Calculating replication for 3 random effects
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL")
{
	rep1<-as.numeric(format(round(repz1/repz2, 0), nsmall=0, scientific=FALSE))
	rep2<-as.numeric(format(round(repz2/repz3, 0), nsmall=0, scientific=FALSE))
	rep0<-repz3
	reps2<-c(rep1, rep2, rep0)
	reps<-reps2
}

#Calculating replication for 2 random effects
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL")
{
	rep1<-as.numeric(format(round(repz1/repz2, 0), nsmall=0, scientific=FALSE))
	rep0<-repz2
	reps1<-c(rep1, rep0)
	reps<-reps1
}




# Calculating the Actual EMS

# One defined random effects
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL")
{
	EMS<- cov12 + rep0*cov11
#	print(EMS)
}

# Two defined random effects
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL")
{
	EMS<-cov23 + rep0*cov22 + rep0*rep2*cov21
#	print(EMS)
}

# Three defined random effects
if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL")
{
	EMS<- cov34 + rep0*cov33 + rep0*rep3*cov32 + rep0*rep3*rep2*cov31
#	print(EMS)
}



#Replication table

HTMLbr()
title<-c("Table of average replication in the original design")
HTML.title(title, HR=2, align="left")
Comment<-c("The following replication of the random effects are used in the power analyses, unless alternative replications are defined by the user. They are an estimate of the replication used within the original design.")
HTML.title("</bf> ", HR=2, align="left")
HTML.title(Comment, HR=0, align="left")

table<-matrix(nrow=1, ncol=length(reffects))
for ( i in 1:length(reffects))
{
	table[1,i]=format(round(reps[i], 0), nsmall=0, scientific=FALSE)
}
vars<-data.frame(table)
colnames(vars)<-reffects
rownames(vars)<-c("Factor replication")
HTML(vars, classfirstline="second", align="left")





#EMS
EMSz<-format(round(EMS, 2), nsmall=2, scientific=FALSE)

HTMLbr()
title<-c("Estimated between subject variance")
HTML.title(title, HR=2, align="left")

Comment<-c("The estimate of the between subject variance is ")
Comment<-paste(Comment, EMSz, sep="")
Comment<-paste(Comment, ". This is computed using the variance components, calculated from the data generated by the original design, and the average replication of the original design.", sep="")
HTML.title("</bf> ", HR=2, align="left")
HTML.title(Comment, HR=0, align="left")



#First plot


HTMLbr()
title<-c("Power curve for an approximate original design")
HTML.title(title, HR=2, align="left")
#print(repz1)
EMS2<-sqrt(EMS)
diffmin<-as.numeric(power.t.test(n=repz1, sd=EMS2, sig.level=siglevel, power=0.05, delta=NULL)[2])
diffmax<-as.numeric(power.t.test(n=repz1, sd=EMS2, sig.level=siglevel, power=0.99, delta=NULL)[2])

diffgap<-(diffmax-diffmin)/10
diffs<-c(diffmin,diffmin+diffgap,diffmin+2*diffgap,diffmin+3*diffgap,diffmin+4*diffgap,diffmin+5*diffgap,diffmin+6*diffgap,diffmin+7*diffgap,diffmin+8*diffgap,diffmin+9*diffgap,diffmin+10*diffgap)




#sampleSize<-c(repz1)
#factorSize<-rep1


# Power analysis functions

# actual change using imputted values
powercurvesactual<-function(standev, sampleSize, factorSize)
{
	legtitle<-c("Replication of subjects = ")
	legtitle2<-rep(legtitle,length(factorSize))
	legtitle3<-paste(legtitle2, factorSize)
	sample<-sampleSize
	temp1<-diffs
	temp2<-matrix(1,nrow=length(diffs),ncol=length(sample))
	for (j in 1:length(sample))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=sample[j], delta=diffs[i], sd=standev, sig.level=siglevel)
	 		temp2[i,j]=test
		}
	}
	temp3<-100*temp2
	powergraph=cbind(diffs,temp3)
	powerPlot <- sub(".html", "powerPlot.jpg", htmlFile)
	jpeg(powerPlot)
	matplot(diffs,temp3, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(sample)), main=graphTitle0, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(sample)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", powerPlot), Align="left")
}



powercurvesactual(EMS, repz1, rep1)


#Replication table

#HTMLbr()
#title<-c("Replication used for approximate original design")
#HTML.title(title, HR=2, align="left")

#table<-matrix(nrow=1, ncol=length(reffects))
#for ( i in 1:length(reffects))
#{
#	table[1,i]=format(round(reps[i], 0), nsmall=0, scientific=FALSE)
#}
#vars<-data.frame(table)
#colnames(vars)<-reffects
#rownames(vars)<-c("Factor replication")
#HTML(vars, classfirstline="second", align="left")




#user defined replication 


if (repList1!="NULL")
{
	tempChanges <-strsplit(repList1, ",")
	userlist1 <- numeric(0)
	for(i in 1:length(tempChanges[[1]])) { userlist1 [length(userlist1)+1] = as.numeric(tempChanges[[1]][i]) } 
#	print(userlist1)
}
if (repList2!="NULL")
{
	tempChanges <-strsplit(repList2, ",")
	userlist2 <- numeric(0)
	for(i in 1:length(tempChanges[[1]])) { userlist2 [length(userlist2)+1] = as.numeric(tempChanges[[1]][i]) } 
#	print(userlist2)
}
if (repList3!="NULL")
{
	tempChanges <-strsplit(repList3, ",")
	userlist3 <- numeric(0)
	for(i in 1:length(tempChanges[[1]])) { userlist3 [length(userlist3)+1] = as.numeric(tempChanges[[1]][i]) } 
#	print(userlist3)
}
if (repList4!="NULL")
{
	tempChanges <-strsplit(repList4, ",")
	userlist4 <- numeric(0)
	for(i in 1:length(tempChanges[[1]])) { userlist4 [length(userlist4)+1] = as.numeric(tempChanges[[1]][i]) } 
#	print(userlist4)
}




# Calculating the EMS's



# Power analysis functions
# actual change using imputted values

#powercurvespredicted<-function(standev, samplesize, userlistz, rdmeffect)
#{
#	legtitle<-paste("- replication of",rdmeffect)
#	legtitle2<-rep(legtitle,length(userlistz))
#	legtitle3<-paste(userlistz,legtitle2)
#	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesize))
#	for (j in 1:length(samplesize))
#	{
#		for(i in 1:length(diffs))
#		{
#			test<-sipowerttest(n=samplesize[j], delta=diffs[i], sd=standev[j], sig.level=siglevel)
#	 		temp[i,j]=test
#		}
#	}
#	temp2<-100*temp
#	powergraph=cbind(diffs,temp2)
#	plot2 <- sub(".html", "plot.jpg", htmlFile)
#	jpeg(plot2)
#	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesize)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
#	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesize)))
#	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot2), Align="left")
#}










# Two random effects

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL" && repList1!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect1, sep="")
	HTML.title(title, HR=2, align="left")

	EMSlist11<-numeric(0)
	samplesizelist11<-numeric(0)
	for(i in 1:length(userlist1))
	{
		EMSlist11[i]<-cov12+ rep0*cov11
		samplesizelist11[i]<-rep0*userlist1[i]
	}



	legtitle<-paste("Replication of ",randomEffect1, "= ")
	legtitle2<-rep(legtitle,length(userlist1))
	legtitle3<-paste(legtitle2, userlist1)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist11))
	for (j in 1:length(samplesizelist11))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist11[j], delta=diffs[i], sd=EMSlist11[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot1 <- sub(".html", "plot1.jpg", htmlFile)
	jpeg(plot1)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist11)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist11)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot1), Align="left")

	Comment1<-c("These power curves assume the replication of the random effect ")
	Comment1<-paste(Comment1, randomEffect2, sep="") 
	Comment1<-paste(Comment1, " is ", sep="") 
	Comment1<-paste(Comment1, rep0, sep="") 
	Comment1<-paste(Comment1, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment1, HR=0, align="left")
}

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL" && repList2!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect2, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect1, sep="")

	HTML.title(title, HR=2, align="left")

	EMSlist12<-numeric(0)
	samplesizelist12<-numeric(0)
	for(i in 1:length(userlist2))
	{
		EMSlist12[i]<-cov12 + userlist2[i]*cov11
		samplesizelist12[i]<-rep1*userlist2[i]
	}


	legtitle<-paste("Replication of ",randomEffect2 , "= ")
	legtitle2<-rep(legtitle,length(userlist2))
	legtitle3<-paste(legtitle2, userlist2)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist12))
	for (j in 1:length(samplesizelist12))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist12[j], delta=diffs[i], sd=EMSlist12[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot2 <- sub(".html", "plot2.jpg", htmlFile)
	jpeg(plot2)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist12)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist12)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot2), Align="left")


	Comment2<-c("These power curves assume the replication of the random effect ")
	Comment2<-paste(Comment2, randomEffect1, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep1, sep="") 
	Comment2<-paste(Comment2, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment2, HR=0, align="left")
}








if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 == "NULL" && randomEffect4 == "NULL" && repList1!="NULL" && repList2!="NULL")
{
	l1<- length(userlist1) 
	l2<- length(userlist2)

	if (l1 ==l2)	
	{
		HTMLbr()
		title<-c("Power curves for a design varying the replication of ")
		title<-paste(title, randomEffect1, sep="")
		title<-paste(title, " and ", sep="")
		title<-paste(title, randomEffect2, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect1, sep="")
		HTML.title(title, HR=2, align="left")
	
		EMSlist12<-numeric(0)
		samplesizelist12<-numeric(0)
		for(i in 1:length(userlist2))
		{
			EMSlist12[i]<-cov12 + userlist2[i]*cov11
			samplesizelist12[i]<-userlist1[i]*userlist2[i]
		}
	
	
	
		legtitle<-paste("Replication of",randomEffect1, "- ")
		legtitle2<-rep(legtitle,length(userlist1))
		legtitle3<-paste(", replication of ",randomEffect2, "- ")
		legtitle4<-rep(legtitle3,length(userlist1))
		legtitle3<-paste(legtitle2, userlist1,  legtitle4, userlist2)


		temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist12))
		for (j in 1:length(samplesizelist12))
		{
			for(i in 1:length(diffs))
			{
				test<-sipowerttest(n=samplesizelist12[j], delta=diffs[i], sd=EMSlist12[j], sig.level=siglevel)
		 		temp[i,j]=test
			}
		}
		temp2<-100*temp
		powergraph=cbind(diffs,temp2)
		plot12 <- sub(".html", "plot12.jpg", htmlFile)
		jpeg(plot12)
		matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist12)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
		legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist12)))
		void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot12), Align="left")
	}
}













# Three random effects

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL" && repList1!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect1)
	HTML.title(title, HR=2, align="left")

	EMSlist31<-numeric(0)
	samplesizelist31<-numeric(0)
	for(i in 1:length(userlist1))
	{
		EMSlist31[i]<-cov23+ rep0*cov22 + rep0*rep2*cov21
		samplesizelist31[i]<-rep0*rep2*userlist1[i]
	}


	legtitle<-paste("Replication of ",randomEffect1, "= ")
	legtitle2<-rep(legtitle,length(userlist1))
	legtitle3<-paste(legtitle2, userlist1)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist31))
	for (j in 1:length(samplesizelist31))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist31[j], delta=diffs[i], sd=EMSlist31[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot1 <- sub(".html", "plot1.jpg", htmlFile)
	jpeg(plot1)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist31)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist31)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot1), Align="left")

	Comment1<-c("These power curves assume the replication of the random effect ")
	Comment1<-paste(Comment1, randomEffect2, sep="") 
	Comment1<-paste(Comment1, " is ", sep="") 
	Comment1<-paste(Comment1, rep2, sep="") 
	Comment1<-paste(Comment1, " and the replication of the random effect ", sep="") 
	Comment1<-paste(Comment1, randomEffect3, sep="") 
	Comment1<-paste(Comment1, " is ", sep="") 
	Comment1<-paste(Comment1, rep0, sep="") 
	Comment1<-paste(Comment1, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment1, HR=0, align="left")
}

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL" && repList2!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect2, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect1, sep="")

	HTML.title(title, HR=2, align="left")

	EMSlist32<-numeric(0)
	samplesizelist32<-numeric(0)
	for(i in 1:length(userlist2))
	{
		EMSlist32[i]<-cov23+ rep0*cov22 + rep0*userlist2[i]*cov21
		samplesizelist32[i]<-rep0*rep1*userlist2[i]
	}


	legtitle<-paste("Replication of ",randomEffect2, "= ")
	legtitle2<-rep(legtitle,length(userlist2))
	legtitle3<-paste(legtitle2, userlist2)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist32))
	for (j in 1:length(samplesizelist32))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist32[j], delta=diffs[i], sd=EMSlist32[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot2 <- sub(".html", "plot2.jpg", htmlFile)
	jpeg(plot2)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist32)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist32)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot2), Align="left")


	Comment2<-c("These power curves assume the replication of the random effect ")
	Comment2<-paste(Comment2, randomEffect1, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep1, sep="") 
	Comment2<-paste(Comment2, " and the replication of ", sep="") 
	Comment2<-paste(Comment2, randomEffect3, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep0, sep="") 
	Comment2<-paste(Comment2, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment2, HR=0, align="left")
}


if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL" && repList3!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect3, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect2, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect1, sep="")
	HTML.title(title, HR=2, align="left")

	EMSlist33<-numeric(0)
	samplesizelist33<-numeric(0)
	for(i in 1:length(userlist3))
	{
		EMSlist33[i]<-cov23+ userlist3[i]*cov22 + userlist3[i]*rep2*cov21
		samplesizelist33[i]<-rep2*rep1*userlist3[i]
	}


	legtitle<-paste("Replication of ",randomEffect3, "= ")
	legtitle2<-rep(legtitle,length(userlist3))
	legtitle3<-paste(legtitle2, userlist3)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist33))
	for (j in 1:length(samplesizelist33))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist33[j], delta=diffs[i], sd=EMSlist33[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot3 <- sub(".html", "plot3.jpg", htmlFile)
	jpeg(plot3)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist33)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist33)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot3), Align="left")


	Comment2<-c("These power curves assume the replication of the random effect ")
	Comment2<-paste(Comment2, randomEffect1, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep1, sep="") 
	Comment2<-paste(Comment2, " and the replication of ", sep="") 
	Comment2<-paste(Comment2, randomEffect2, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep2, sep="") 
	Comment2<-paste(Comment2, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment2, HR=0, align="left")
}











if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 == "NULL" && repList1!="NULL" && repList2!="NULL"&& repList3!="NULL")
{
	l1<- length(userlist1) 
	l2<- length(userlist2)
	l3<- length(userlist3)

	if (l1 ==l2 && l1==l3)	
	{
		HTMLbr()
		title<-c("Power curves for a design varying the replication of ")
		title<-paste(title, randomEffect1, sep="")
		title<-paste(title, ", ", sep="")
		title<-paste(title, randomEffect2, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect1, sep="")

		title<-paste(title, " and ", sep="")
		title<-paste(title, randomEffect3, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect2, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect1, sep="")

		HTML.title(title, HR=2, align="left")
	
		EMSlist12<-numeric(0)
		samplesizelist12<-numeric(0)
		for(i in 1:length(userlist2))
		{
			EMSlist12[i]<-cov23+ userlist3[i]*cov22 + userlist3[i]*userlist2[i]*cov21
			samplesizelist12[i]<-userlist1[i]*userlist2[i]*userlist3[i]
		}
	
	
	
		legtitle<-paste("Replication of",randomEffect1, "- ")
		legtitle2<-rep(legtitle,length(userlist1))
		legtitle3<-paste(", replication of ",randomEffect2, "- ")
		legtitle4<-rep(legtitle3,length(userlist1))
		legtitle5<-paste(", replication of ",randomEffect3, "- ")
		legtitle6<-rep(legtitle5,length(userlist1))
		legtitle3<-paste(legtitle2, userlist1,  legtitle4, userlist2,  legtitle6, userlist3)


		temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist12))
		for (j in 1:length(samplesizelist12))
		{
			for(i in 1:length(diffs))
			{
				test<-sipowerttest(n=samplesizelist12[j], delta=diffs[i], sd=EMSlist12[j], sig.level=siglevel)
		 		temp[i,j]=test
			}
		}
		temp2<-100*temp
		powergraph=cbind(diffs,temp2)
		plot12 <- sub(".html", "plot12.jpg", htmlFile)
		jpeg(plot12)
		matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist12)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
		legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist12)))
		void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot12), Align="left")
	}
}





























# Three random effects

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL" && repList1!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect1, sep="")
	HTML.title(title, HR=2, align="left")

	EMSlist41<-numeric(0)
	samplesizelist41<-numeric(0)
	for(i in 1:length(userlist1))
	{
		EMSlist41[i]<-cov34 + rep0*cov33 + rep0*rep3*cov32 + rep0*rep3*rep2*cov31
		samplesizelist41[i]<-rep0*rep2*rep3*userlist1[i]
	}

	legtitle<-paste("Replication of ",randomEffect1, "= ")
	legtitle2<-rep(legtitle,length(userlist1))
	legtitle3<-paste(legtitle2,userlist1)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist41))
	for (j in 1:length(samplesizelist41))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist41[j], delta=diffs[i], sd=EMSlist41[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot1 <- sub(".html", "plot1.jpg", htmlFile)
	jpeg(plot1)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist41)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist41)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot1), Align="left")

	Comment1<-c("These power curves assume the replication of the random effect ")
	Comment1<-paste(Comment1, randomEffect2, sep="") 
	Comment1<-paste(Comment1, " is ", sep="") 
	Comment1<-paste(Comment1, rep2, sep="") 
	Comment1<-paste(Comment1, ", the replication of the random effect ", sep="") 
	Comment1<-paste(Comment1, randomEffect3, sep="") 
	Comment1<-paste(Comment1, " is ", sep="") 
	Comment1<-paste(Comment1, rep3, sep="") 
	Comment1<-paste(Comment1, " and the replication of the random effect ", sep="") 
	Comment1<-paste(Comment1, randomEffect4, sep="") 
	Comment1<-paste(Comment1, " is ", sep="") 
	Comment1<-paste(Comment1, rep0, sep="") 
	Comment1<-paste(Comment1, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment1, HR=0, align="left")
}

if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL" && repList2!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect2, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect1, sep="")

	HTML.title(title, HR=2, align="left")

	EMSlist42<-numeric(0)
	samplesizelist42<-numeric(0)
	for(i in 1:length(userlist2))
	{
		EMSlist42[i]<-cov34 + rep0*cov33 + rep0*rep3*cov32 + rep0*rep3*userlist2[i]*cov31
		samplesizelist42[i]<-rep0*rep1*rep3*userlist2[i]
	}

	legtitle<-paste("Replication of ",randomEffect2 , "= ")
	legtitle2<-rep(legtitle,length(userlist2))
	legtitle3<-paste(legtitle2,userlist2)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist42))
	for (j in 1:length(samplesizelist42))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist42[j], delta=diffs[i], sd=EMSlist42[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot2 <- sub(".html", "plot2.jpg", htmlFile)
	jpeg(plot2)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist42)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist42)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot2), Align="left")


	Comment2<-c("These power curves assume the replication of the random effect ")
	Comment2<-paste(Comment2, randomEffect1, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep1, sep="") 
	Comment2<-paste(Comment2, ", the replication of the random effect ", sep="") 
	Comment2<-paste(Comment2, randomEffect3, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep3, sep="") 
	Comment2<-paste(Comment2, " and the replication of ", sep="") 
	Comment2<-paste(Comment2, randomEffect4, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep0, sep="") 
	Comment2<-paste(Comment2, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment2, HR=0, align="left")
}


if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL" && repList3!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect3, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect2, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect1, sep="")
	HTML.title(title, HR=2, align="left")

	EMSlist43<-numeric(0)
	samplesizelist43<-numeric(0)
	for(i in 1:length(userlist3))
	{
		EMSlist43[i]<-cov34 + rep0*cov33 + rep0*userlist3[i]*cov32 + rep0*userlist3[i]*rep2*cov31
		samplesizelist43[i]<-rep2*rep1*rep0*userlist3[i]
	}


	legtitle<-paste("Replication of ",randomEffect3, "= ")
	legtitle2<-rep(legtitle,length(userlist3))
	legtitle3<-paste(legtitle2,userlist3)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist43))
	for (j in 1:length(samplesizelist43))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist43[j], delta=diffs[i], sd=EMSlist43[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot3 <- sub(".html", "plot3.jpg", htmlFile)
	jpeg(plot3)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist43)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist43)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot3), Align="left")


	Comment2<-c("These power curves assume the replication of the random effect ")
	Comment2<-paste(Comment2, randomEffect1, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep1, sep="") 
	Comment2<-paste(Comment2, ", the replication of the random effect ", sep="") 
	Comment2<-paste(Comment2, randomEffect2, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep2, sep="") 
	Comment2<-paste(Comment2, " and the replication of ", sep="") 
	Comment2<-paste(Comment2, randomEffect4, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep0, sep="") 
	Comment2<-paste(Comment2, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment2, HR=0, align="left")
}





if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL" && repList4!="NULL")
{
	HTMLbr()
	title<-c("Power curves for a design varying the replication of ")
	title<-paste(title, randomEffect4, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect3, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect2, sep="")
	title<-paste(title, " within ", sep="")
	title<-paste(title, randomEffect1, sep="")
	HTML.title(title, HR=2, align="left")

	EMSlist44<-numeric(0)
	samplesizelist44<-numeric(0)
	for(i in 1:length(userlist4))
	{
		EMSlist44[i]<-cov34 + rep0*cov33 + userlist4[i]*rep3*cov32 + rep3*userlist4[i]*rep2*cov31
		samplesizelist44[i]<-rep2*rep1*rep3*userlist4[i]
	}


	legtitle<-paste("Replication of ",randomEffect4, "= ")
	legtitle2<-rep(legtitle,length(userlist4))
	legtitle3<-paste(legtitle2,userlist4)
	temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist44))
	for (j in 1:length(samplesizelist44))
	{
		for(i in 1:length(diffs))
		{
			test<-sipowerttest(n=samplesizelist44[j], delta=diffs[i], sd=EMSlist44[j], sig.level=siglevel)
	 		temp[i,j]=test
		}
	}
	temp2<-100*temp
	powergraph=cbind(diffs,temp2)
	plot4 <- sub(".html", "plot4.jpg", htmlFile)
	jpeg(plot4)
	matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist44)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
	legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist44)))
	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot4), Align="left")


	Comment2<-c("These power curves assume the replication of the random effect ")
	Comment2<-paste(Comment2, randomEffect1, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep1, sep="") 
	Comment2<-paste(Comment2, ", the replication of the random effect ", sep="") 
	Comment2<-paste(Comment2, randomEffect2, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep2, sep="") 
	Comment2<-paste(Comment2, " and the replication of ", sep="") 
	Comment2<-paste(Comment2, randomEffect3, sep="") 
	Comment2<-paste(Comment2, " is ", sep="") 
	Comment2<-paste(Comment2, rep3, sep="") 
	Comment2<-paste(Comment2, ".", sep="") 
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(Comment2, HR=0, align="left")
}









if (randomEffect1 != "NULL" && randomEffect2 != "NULL" && randomEffect3 != "NULL" && randomEffect4 != "NULL" && repList1!="NULL" && repList2!="NULL" && repList3!="NULL" && repList4!="NULL")
{
	l1<- length(userlist1) 
	l2<- length(userlist2)
	l3<- length(userlist3)
	l4<- length(userlist4)

	if (l1 ==l2 && l1==l3 && l1==l4)	
	{
		HTMLbr()
		title<-c("Power curves for a design varying the replication of ")
		title<-paste(title, randomEffect1, sep="")
		title<-paste(title, ", ", sep="")
		title<-paste(title, randomEffect2, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect1, sep="")

		title<-paste(title, ", ", sep="")
		title<-paste(title, randomEffect3, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect2, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect1, sep="")

		title<-paste(title, " and ", sep="")
		title<-paste(title, randomEffect4, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect3, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect2, sep="")
		title<-paste(title, " within ", sep="")
		title<-paste(title, randomEffect1, sep="")
		HTML.title(title, HR=2, align="left")
	
		EMSlist12<-numeric(0)
		samplesizelist12<-numeric(0)
		for(i in 1:length(userlist2))
		{
			EMSlist12[i]<-cov34 + userlist4[i]*cov33 + userlist4[i]*userlist3[i]*cov32 + userlist4[i]*userlist3[i]*userlist2[i]*cov31
			samplesizelist12[i]<-userlist1[i]*userlist2[i]*userlist3[i]*userlist4[i]
		}
	
	
	
		legtitle<-paste("Rep. of",randomEffect1, "- ")
		legtitle2<-rep(legtitle,length(userlist1))
		legtitle3<-paste(", rep. of ",randomEffect2, "- ")
		legtitle4<-rep(legtitle3,length(userlist1))
		legtitle5<-paste(", rep. of ",randomEffect3, "- ")
		legtitle6<-rep(legtitle5,length(userlist1))
		legtitle7<-paste(", rep. of ",randomEffect4, "- ")
		legtitle8<-rep(legtitle5,length(userlist1))
		legtitle3<-paste(legtitle2, userlist1,  legtitle4, userlist2,  legtitle6, userlist3,  legtitle8, userlist4)


		temp<-matrix(1,nrow=length(diffs),ncol=length(samplesizelist12))
		for (j in 1:length(samplesizelist12))
		{
			for(i in 1:length(diffs))
			{
				test<-sipowerttest(n=samplesizelist12[j], delta=diffs[i], sd=EMSlist12[j], sig.level=siglevel)
		 		temp[i,j]=test
			}
		}
		temp2<-100*temp
		powergraph=cbind(diffs,temp2)
		plot12 <- sub(".html", "plot12.jpg", htmlFile)
		jpeg(plot12)
		matplot(diffs,temp2, lwd=2,pch = 1: length(temp2), ylim=c(0, 100), type = "l",  col = rainbow(length(samplesizelist12)), main=graphTitle1, xlab = "Biological difference", ylab="Power (%)")
		legend("bottomright",legend= legtitle3, cex=0.6,lwd=2,lty=1:length(temp2),bg="white",col=rainbow(length(samplesizelist12)))
		void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plot12), Align="left")
	}
}



#References

HTMLbr()
HTML.title("<bf>Statistical references", HR=2, align="left")


HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>  Snedecor GW and Cochran WG. (1989). Statistical Methods (8th Edition). Iowa State University Press, Ames, Iowa.", HR=0, align="left")


HTMLbr()
HTML.title("<bf>R references", HR=2, align="left")
	
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   R Development Core Team (2008). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.", HR=0, align="left")

#nlme
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> 
	Jose Pinheiro, Douglas Bates, Saikat DebRoy, Deepayan Sarkar and the R Core team (2008). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-90.
	", HR=0, align="left")

#R2HTML
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
	Lecoutre, Eric (2003). The R2HTML Package. R News, Vol 3. N. 3, Vienna, Austria.
	", HR=0, align="left")
